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A  two-dimensional  two-phase  mass  transport  model  for  liquid-feed  direct  methanol  fuel  cells  (DMFCs) 
is  presented  in  this  paper.  The  fluid  flow  and  mass  transport  across  the  membrane  electrode  assembly 
(MEA)  is  formulated  based  on  the  classical  multiphase  flow  theory  in  the  porous  media.  The  modeling  of 
mass  transport  in  the  catalyst  layers  (CLs)  and  membrane  is  given  more  attentions.  The  effect  of  the  two- 
dimensional  migration  of  protons  in  the  electrolyte  phase  on  the  liquid  flow  behavior  is  considered.  Water 
and  methanol  crossovers  through  the  membrane  are  implicitly  calculated  in  the  governing  equations  of 
momentum  and  methanol  concentration.  A  modified  agglomerate  model  is  developed  to  characterize  the 
microstructure  of  the  CLs.  A  self-written  computer  code  is  used  to  solve  the  inherently  coupled  differential 
governing  equations.  Then  this  model  is  applied  to  investigate  the  mechanisms  of  species  transport  and 
the  distributions  of  the  species  concentrations,  overpotential  and  the  electrochemical  reaction  rates  in 
CLs.  The  effects  of  radius  and  overlapping  angle  of  agglomerates  on  cell  performance  are  also  explored  in 
this  work. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  liquid-feed  direct  methanol  fuel  cell  (DMFC)  is  presently 
regarded  as  one  of  the  most  promising  candidates  to  replace  the 
conventional  batteries,  especially  those  used  for  portable  elec¬ 
tronic  devices  [1].  However,  the  commercialization  of  DMFCs  is 
still  hindered  by  several  technological  problems  [2],  such  as  low 
electro-activity  of  methanol  oxidation  on  the  surface  of  anode  cat¬ 
alyst  particles,  substantial  methanol  crossover  through  the  proton 
exchange  membrane  (MEM)  from  the  anode  to  the  cathode  and  the 
water  management  to  prevent  severe  cathode  flooding.  Extensive 
efforts  have  been  made  to  solve  these  problems  and  also  much  work 
has  been  focused  on  the  fundamentals  of  the  DMFC  system,  among 
which  the  computational  modeling  of  DMFCs  played  an  important 
role.  Modeling  was  regarded  as  a  powerful  and  economical  tool  to 
investigate  the  intrinsically  coupled  processes  in  a  DMFC  which  are 
difficult  to  be  studied  by  experiments. 

Over  the  past  decade,  many  mathematical  models  have  been 
developed  [3-19],  which  include  single-phase  models  [3-11]  and 
two-phase  models  [12-19].  Most  of  the  previous  published  mod¬ 
els  for  DMFCs  are  based  on  the  single-phase  flow  assumption  and 
several  significant  conclusions  have  been  obtained  according  to 
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the  analysis  of  the  distribution  of  species  concentration,  the  effect 
of  operation  conditions  on  cell  performance  and  so  on.  However, 
the  in  situ  visualization  research  of  DMFC  [20-23]  revealed  that 
the  fluid  flow  in  the  DMFC  is  not  a  simple  single-phase  process, 
there  is  coexisting  liquid  and  gas  two-phase  flow  behavior  in  the 
flow  channel  because  of  the  generation  of  gaseous  carbon  dioxide 
in  the  anode  catalyst  layer  (ACL)  and  formation  of  liquid  water  in 
the  cathode  catalyst  layer  (CCL).  The  significant  effects  of  the  two- 
phase  flow  behavior  on  the  mass  transport  processes  and  on  the 
cell  performance  [12,17,19]  demand  the  development  of  two-phase 
mathematical  models. 

Most  of  the  published  two-phase  models  are  based  on  the  clas¬ 
sical  multiphase  flow  theory  with  assumptions  and  simplifications. 
Wang  and  Wang  [12]  developed  a  two-phase  model  for  DMFCs 
based  on  the  mixture  multiphase  flow  approach  for  the  porous 
regions  and  the  drift  flux  model  for  the  two-phase  flow  in  channels 
of  the  flow-field  plates.  The  liquid  and  gas  mixture  is  assumed  to  be 
a  homogeneous  fluid  at  the  thermodynamic  equilibrium  condition, 
and  its  properties  can  be  derived  from  the  parameters  of  each  indi¬ 
vidual  phase  based  on  certain  mixing  rules.  The  CLs  in  this  model 
are  treated  as  infinitely  thin  interfaces.  Ge  and  Liu  [13]  presented  a 
three-dimensional,  two-phase  mathematical  model  adopting  the 
general  momentum  conservation  equations  with  a  source  term 
accounting  for  the  forces  exerted  on  the  fluid  by  the  solid  matrix  of 
the  porous  media.  The  CLs  are  treated  as  porous  regions  with  finite 
width.  Differential  equations  are  solved  by  the  CFD  technique  on 
the  single  computational  domain  consisting  of  seven  layers:  flow 
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channels,  diffusion  layers  (DLs),  CLs  and  the  MEM.  Divisek  et  al. 
[  14]  developed  a  two-dimensional,  two-phase  mathematical  model 
for  a  DMFC  by  regarding  the  DLs  as  water-gas  systems  in  the  pore 
space  with  saturation  permeability  varying  according  to  the  cap¬ 
illary  effects.  A  new  parameterization  of  the  relationship  between 
capillary  pressure  and  saturation  is  introduced  in  their  work.  In 
particular,  the  electrochemical  reactions  are  split  up  into  reaction 
chains  to  obtain  a  more  precise  expression  of  the  electrochemical 
reaction  rates. 

From  the  models  mentioned  above,  it  can  be  concluded  that 
the  mass  transport  processes  in  the  CLs  and  in  the  MEM  are  not 
fully  explored.  Many  researchers  treated  the  CLs  as  either  infinitely 
thin  interfaces  or  pseudo-homogeneous  porous  media  to  simplify 
their  models.  The  microstructure  of  the  CLs  is  ignored  and  the 
simulations  of  potential  and  velocities  in  electrolyte-phase-region 
(CLs  and  MEM)  are  simplified.  Studies  of  the  microstructure  of 
CLs  by  high-resolution  scanning  electron  microscopy  (HR-SEM) 
[24]  showed  that  the  CLs  consist  of  randomly  distributed  pores 
and  agglomerates  of  much  smaller  carbon-based  catalyst  particles. 
Computational  model  invoking  agglomerate  approach  to  model  the 
mass  transport  and  electrochemical  reaction  processes  in  the  CLs 
is  first  developed  for  the  alkaline  oxygen  electrode  by  Giner  and 
Hunter  [25]  with  the  concept  of  “flooded  agglomerates”.  In  the  case 
of  PEM  fuel  cells,  researchers  have  also  studied  various  phenomena 
in  the  CLs  based  on  the  agglomerate  model  [26-30]. 

However,  only  a  few  literatures  adopted  the  agglomerate  model 
to  characterize  the  CLs  of  the  DMFCs.  Nordlund  and  Lindbergh  [31  ] 
modeled  the  anode  of  a  DMFC  based  on  an  agglomerate  model 
for  the  mass  transport  and  electrochemical  reaction  processes  of 
methanol  and  carbon  dioxide  in  the  ACL.  Modeling  results  showed 
that  the  concentration  profiles  within  the  agglomerates  on  the 
anode  can  be  neglected  so  that  the  model  can  be  simplified.  In  their 
model,  they  assumed  isobaric  condition  in  the  electrode  and  that 
carbon  dioxide  is  dissolved  in  the  water.  Yang  and  Zhao  [15]  pre¬ 
sented  a  two-dimensional,  two-phase  model  for  liquid-feed  DMFCs 
with  a  micro-agglomerate  model  developed  for  the  CCL.  The  oxy¬ 
gen  transfer  processes  from  the  gas  pores  to  the  active  sites  were 
modeled  by  considering  the  agglomerate  covered  by  a  liquid  film. 

Some  simplifications  used  in  the  two-phase  mass  transport 
models  above,  including  models  that  adopted  the  agglomerate 
approach,  can  be  collected  as  follows:  (1)  The  overpotential  in 
the  CLs  is  always  assumed  to  distribute  uniform  as  the  proton 
migration  process  in  the  electrolyte  phase  of  CLs  is  not  consid¬ 
ered.  The  electrochemical  reaction  rates  are  very  sensitive  to  the 
overpotential,  so  it  is  worth  investigating  the  effect  of  overpoten¬ 
tial  on  the  distribution  of  electrochemical  reaction  rates  in  the 
CLs;  (2)  Both  water  and  methanol  crossover  rates  are  explicitly 
calculated  by  expressions  that  account  for  the  combination  of  diffu¬ 
sion,  back  convection  and  electro-osmotic  drag.  Consequently,  the 
two-dimensional  dispersion  of  water  and  methanol  in  the  MEM 
is  simplified  as  one-dimensional  permeation  process  through  the 
MEM  because  only  the  variables  at  the  two  boundaries,  interfaces 
between  the  CLs  and  the  MEM,  can  be  used  in  these  expressions;  (3) 
The  ‘parasitic’  methanol  current  density  is  dictated  by  the  methanol 
crossover  rate.  Since  the  mass  transport  of  methanol  in  the  CCL  is 
not  considered,  the  distribution  of  the  ‘parasitic’  current  cannot 
be  obtained;  (4)  The  geometry  of  the  spherical  agglomerate  with  a 
Nafion  or  water  coating  is  always  adopted  in  the  agglomerate  mod¬ 
els.  However,  the  Nafion  coating  or  water  film  covering  outside  of 
the  spherical  agglomerate  will  block  the  access  of  electrons,  so  a 
more  realistic  geometry  of  the  agglomerate  should  be  developed  to 
satisfy  the  common  recognition  that  the  CLs  should  supply  access 
for  all  the  three  phases:  reactants,  protons  and  electrons. 

The  objective  of  this  work  is  to  develop  a  two-dimensional, 
two-phase  mass  transport  model  for  a  liquid  feed  DMFC  which 


can  provide  a  useful  insight  into  the  mechanisms  of  species  trans¬ 
port  and  more  realistic  distributions  of  variable  fields  in  the  CLs. 
In  this  model,  the  local  overpotential  in  the  CLs  are  derived  from 
the  current  conservation  equations  by  considering  the  transport  of 
electrons  in  the  carbon  phase  and  protons  in  the  electrolyte  phase. 
Water  and  methanol  crossovers  through  the  MEM  are  taken  into 
account  by  incorporating  the  MEM  into  the  computational  domains 
of  the  momentum  conservation  equations  and  the  species  equation 
for  methanol.  The  ‘parasitic’  current  density  in  the  CCL  is  calculated 
by  a  Tafel-like  expression  similar  to  that  in  the  anode  side  and  a 
modified  agglomerate  model  is  applied  for  the  CCL. 

2.  Mathematical  model 

The  two-dimensional  computational  domain,  as  enclosed  by  the 
dashed  lines  in  Fig.  1,  represents  the  geometry  of  a  complete  MEA, 
which  consists  of  anode  diffusion  layer  (ADL),  ACL,  MEM,  CCL  and 
cathode  diffusion  layer  (CDL).  Since  the  parallel  flow-field  plates  are 
used  in  the  DMFC,  only  a  periodic  unit:  part  of  the  MEA  sandwiched 
by  a  flow  channel  and  two  hemi-ribs  is  considered  in  this  model. 
Governing  equations  are  elaborated  below. 

2  A.  Potential  governing  equations 

In  a  DMFC,  protons  and  electrons  are  generated  in  the  ACL  and 
combined  in  the  CCL.  The  motion  of  protons  and  electrons  obeys 
the  Ohm’s  law  and  the  current  conservation  equation: 

1  =  V(cr<p)  (1) 

V ■ I = i  (2) 

where  1  is  the  current  vector  and  i  is  the  generation  rate  of  cur¬ 
rent;  a  is  the  electric  conductivity  of  materials  and  cp  is  the  electric 
potential. 

Following  the  concept  of  electric  potential  of  the  membrane 
phase  (pm  (or  electrolyte  phase)  and  solid  phases  <ps  (or  carbon 
phase)  introduced  by  Kulikovs ky  [6],  we  can  obtain  potential  gov¬ 
erning  equations  by  substituting  Eq.  (1)  into  Eq.  (2): 


V  •  (crmv^m)  —  im  (3) 

V -{crsW(ps)  =  is  (4) 

where  crm  and  crs  denote  the  conductivity  of  electrolyte  and  car¬ 
bon  phases;  im  and  is  represent  the  generation  rates  of  protons  and 
electrons  respectively.  In  this  model,  the  ACL,  MEM  and  CCL  are 
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Fig.  1.  Schematic  of  the  computational  domain  of  a  liquid-feed  DMFC. 
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specified  as  the  electrolyte-phase-region,  which  means  the  com¬ 
putational  domain  of  electrolyte  phase  potential. 

2.2.  Velocity  and  pressure  equations  for  gas  and  liquid  phases 

In  the  classical  two-phase  flow  theory  in  a  porous  media,  the 
mass  conservation  equations  can  be  given  as  follows: 


2.3.  Liquid  saturation  in  the  porous  media 

The  liquid  pressure  in  the  porous  media  is  not  equal  to  the  gas 
pressure  because  of  the  capillary  effect.  The  pressure  difference 
between  these  two  phases  is  defined  as  capillary  pressure.  The 
Leverette  function  [32,33]  is  adopted  to  describe  the  relationship 
between  the  liquid  saturation  and  the  capillary  pressure. 


V  ■  (PgUg)  =  rhg  (5) 

v-(piui)  =  mi  (6) 

where  p,  u  and  m  represent  the  density,  the  superficial  velocity 
vector  and  the  mass  generation  rate. 


/  £  \  0*5 

Pc=Pg-P\  =  a  COS  6cl  -)  J(s) 


(  1 .417(1  —  s)  -  2.12(1  +1.263(1  0°  <  0C  <  90° 

\  1.417s  -  2.12s2  +  1.263s3  90°  <  6C  <  180° 


(13) 

(14) 


2.2.1.  Velocity  and  pressure  of  the  gas  phase 

Gas  phase  in  the  DMFC  is  a  multi-component  mixture,  includ¬ 
ing  carbon  dioxide,  water  vapor,  methanol  vapor  and  oxygen  when 
pure  oxygen  is  supplied  in  the  cathode  flow  channel  (CFC).  In  order 
to  simplify  the  model,  the  commonly  used  assumptions  are  also 
adopted  in  the  present  model  that  the  transport  of  methanol  in 
gas  phase  is  neglected  [13,15,18]  and  carbon  dioxide  produced  by 
the  methanol  oxidation  reaction  (MOR)  on  the  anode  and  methanol 
direct  oxidation  reaction  (MDOR)  on  the  cathode  only  exists  in  the 
gas  phase  [15-19].  The  convective  transport  of  gases  is  caused  by 
the  pressure  gradient  of  the  gas  phase,  so  the  expression  of  gas 
velocity  is  given  as  [15]: 

ug  =  -K^Vpg  (7) 

Mg 

where  I<  represents  the  absolute  permeability  of  the  porous  media, 
krg  represents  the  relative  permeability  of  gas  phase,  /xg  and  pg  are 
the  viscosity  and  pressure  of  the  gas  phase.  By  substituting  Eq.  (7) 
into  Eq.  (5),  the  governing  equation  of  gas  phase  pressure  can  be 
obtained  in  the  following  form: 


where  a,  0C  and  £  denote  the  surface  tension  of  the  liquid  phase, 
the  contact  angle  between  the  liquid  phase  and  the  solid  wall,  and 
the  porosity  of  the  porous  media  respectively.  Substituting  Eq.  (13) 
into  Eqs.  (8)  and  (10)  yields: 


Vs 


=  m{  + 


nd^w , 


A^riMg  ^ 
PgkrgMl  S 


(ADL,  ACL) 


(15) 


PgkTgfi\ 

P\K  lMg 


(£H 

(■  ttdMw  .  \ 


(CDL,  CCL) 


(16) 


In  order  to  guarantee  the  robustness  of  the  self-written  code 
and  smaller  interpolation  error,  governing  equations  for  liquid 
saturation  are  restricted  to  the  specified  computational  domains 
indicated  in  Eqs.  (15)  and  (16). 


V  • 


-I< 


Pgkrg 

Mg 


VPg 


=  /T7o 


(8) 


2.2.2.  Velocity  and  pressure  of  the  liquid  phase 

In  a  liquid-feed  DMFC,  the  dilute  methanol  aqueous  solution  is 
supplied  in  the  anode  flow  channel  (AFC).  As  protons  produced  on 
the  anode  exist  in  the  form  of  hydrated  ions,  the  migration  of  pro¬ 
tons  in  an  electrical  field  will  lead  to  the  bulk  motion  of  liquid  phase. 
So  it  can  be  seen  that  not  only  the  gradient  of  liquid  pressure  but 
also  electro-osmotic  drag  force  contributes  to  the  liquid  velocity. 
The  expression  of  liquid  velocity  similar  to  that  used  in  literature 
[12]  is  applied: 


„kr  i  ndM  1 

Uj  =  —K  —  Vpj  h — - —  — 
Ml  A  F 


(9) 


where  nd  represents  the  electro-osmotic  drag  coefficient  of  water.  F 
is  the  Faraday’s  constant.  M  is  equal  to  the  water  molecular  weight 
Mw  due  to  the  dilute  solution  assumption.  Substituting  Eq.  (9)  into 
Eq.  (6),  with  the  aid  of  Eqs.  (l)-(3),  we  obtain: 


V  • 


=  fnx- 


ndMw . 


(10) 


The  relative  permeability  of  gas  phase  krg  and  liquid  phase  krX  in 
Eqs.  (7)-(10)  can  be  formulated  by  the  liquid  saturation  alone  [12]: 


kri=s3  (11) 

krg  =  (1  —  S)3  (12) 


The  mass  conservation  and  momentum  conservation  in  the 
DMFC  can  be  satisfied  by  solving  Eqs.  (7)-(10)  iteratively. 


2.4.  Concentration  of  species 

We  now  turn  our  attention  to  the  mass  transport  of  different 
species,  including  methanol,  oxygen  and  water  vapor,  in  the  porous 
media.  The  general  species  conservation  equations  [34]  can  be  writ¬ 
ten  as  follows: 


V  •  (UjC))  =  V  •  (D?ffVQ)  +  Rj  (17) 

where  Q  and  Rz  represent  the  concentration  and  the  molar  gen¬ 
eration  rate  of  species  i.  Dff  is  the  effective  diffusion  coefficient  of 
species  i  corrected  with  the  porosity  and  liquid  saturation.  uz  stands 
for  the  velocity  of  liquid  or  gas  phase  and  can  be  set  as  the  value 
below  for  different  species: 


/  Ul  CH3OH 

1  \ug  02,H20t 

where  H2Ot  represents  the  water  vapor. 


(18) 


2.5.  Electrochemical  kinetics 


As  the  mechanisms  of  the  complicated  multi-step  electrochem¬ 
ical  reactions  occurring  in  a  DMFC  have  not  been  completely 
understood,  the  Tafel-like  expressions  [18]  are  used  for  the  kinetics 
of  MOR  on  the  anode  and  oxygen  reduction  reaction  (ORR)  on  the 
cathode  for  simplification. 


;  _  ,4  ireff  CM§M 

ia  -  s/taiM  1  -^r 
V 


Y 

exp 


ocaE 

~RT 


(<Fs,a  — 


^m,a) 


(19) 
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Fig.  2.  Sketch  of  the  agglomerate  geometry  in  the  CCL. 


I'c  =  (1  -  s)Aci’o2f 


exp 


~acF 

~RT 


(^m,c 


(20) 


where  CM  and  C02  are  the  concentrations  of  methanol  and  oxygen  in 
the  bulk  pores  of  the  CLs;  C^f  and  C^f  are  reference  concentrations; 
§M  and  §o2  are  correction  factors  derived  from  the  agglomerate 
model  and  will  be  discussed  in  detail  later.  The  reaction  order  y  is 
determined  by  the  methanol  concentration  at  the  active  sites  of  the 
catalyst  particles,  CM§ m,  and  can  be  set  as: 


f  o  cm§m  >  qf 
\  1  Cm£m  <  qf 


(21) 


The  correction  factors,  and  §o2  in  Eqs.  (19)  and  (20),  consider 
effects  of  the  mass  transfer  resistance  induced  by  agglomerates 
with  a  Nation  coating  and  a  overlapping  region  between  two  adja¬ 
cent  agglomerates  as  illustrated  in  Fig.  2.  In  the  ACL,  the  Thiele 
module  is  very  small  so  that  there  is  no  need  to  account  for  the 
concentration  difference  in  the  agglomerate  [31  ].  In  this  model,  the 
correction  factor  £M  =  1.0  is  adopted  for  the  ACL. 

In  the  CCL,  oxygen  in  the  bulk  pores  should  first  dissolves  into 
the  Nation  coating.  This  process  can  be  described  by  the  Henry’s 
law: 


c°.«  -  v  <221 

where  kH  is  the  Henry’s  constant  and  Co2n  is  the  concentration  of 
dissolved  oxygen  at  the  outer  surface  of  the  Nation  coating.  It  is 
assumed  that  the  transport  of  oxygen  in  the  Nation  coating  and  in 
the  agglomerate  is  a  diffusion  dominating  process,  so  that  the  mass 
transport  of  oxygen  through  the  Nation  coating  and  consumption 
of  oxygen  in  the  agglomerate  can  be  modeled  with: 


d2Co2,N  1  dCo2 5n 
dr2  +  r  dr 


(Sub-domain  I) 


(23) 


d2Co2;agg  1^  dCp2 ;  agg 

dr2  r  dr 


—ff — Co2,agg  (Sub-domain  II) 

^02,agg 


(24) 


d^agg  +  2  dC^agg  =  — Co2,agg  (Sub-domain  III)  (25) 

^02,agg 

where  Co2,n  and  C02,agg  represent,  respectively,  the  concentration 
of  oxygen  in  the  Nation  coating  and  in  the  agglomerate.  agg  is 
the  effective  diffusion  coefficient  of  oxygen  in  the  agglomerate  and 
kv  is  the  volumetric  reaction  rate  constant  in  the  agglomerate.  The 
detailed  derivation  of  Eqs.  (23)-(25)  is  enclosed  in  Appendix  A.  The 
correction  factor  §o2  can  be  obtained  by  solving  Eqs.  (23)-(25).  The 
expression  of  the  correction  factor  and  related  correlations  in  detail 
are  listed  in  Table  1. 


2.6.  Current  balance  and  cell  voltage 


In  a  working  fuel  cell,  which  is  part  of  the  circuit,  the  mean 
current  densities  at  the  anode  and  cathode  can  be  respectively 
calculated  by: 


fj^dxdy 

hc  +  hr 


(26) 


Table  1 

Expressions  of  correction  factor  and  related  correlations  for  the  agglomerate  model 


Descriptions 


Expressions 


Total  correction  factor 

Correction  factor  for  Nafion 
coating 

Correction  factor  for 
agglomerate 

The  Thiele  modulus 

The  volume  reaction  rate 
constant 

Function  related  to  overlapping  f(0)  = 
angle 


Function  related  to  the  Thiele 
modulus 


£o2  -  j^N^agg 

£n  = 


l+(2/3)(D^fg/DN)(0s  COth  0s -1)  ln(J?agg+<5N)/Ragg 
6  R2 


i'(4> s) 

:  /(0s 

0S  =  R[n  y/kv/D\ 

/<V  = 


eff 

02,agg 


_ ecc)De»  Cref  “P  [  fa  {<Pm  ~  <Pc )J 


4F(1 


02,agg  02 

cos  9 


[2-(l-cos  0f( 2+cos0)] 


/(0s)  —  ^  JC^(^agg  -  ^in)n] 
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Fig.  3.  Specified  computational  sub-domains  for  the  solving  variables. 


r  ffccL^dxdy 

c  hc  +  h 


Because  part  of  the  oxygen  is  consumed  by  the  MDOR  in  the  CCL, 
the  current  conservation  equation  can  be  expressed  as: 


Alell  —  fa  —  fc  —  fp 


(28) 


where  the  ‘parasitic’  current,  7P,  is  a  virtual  current  and  denotes  the 
fuel  loss  caused  by  methanol  crossover.  It  is  worth  to  point  out  that, 
unlike  the  assumptions  adopted  by  many  other  literatures  that  the 
‘parasitic’  current  is  distributed  uniformly  on  the  x  direction  in  the 
CCL,  the  ‘parasitic’  current  in  this  model  is  calculated  by  adopting 
a  Tafel-like  expression  similar  to  that  used  in  literature  [6]. 


ip  -  Aci$ 


(  Cm£m 

v  qf 


Y 

exp 


Oaf 

~w 


(<^m,c  —  <Ps,c) 


(29) 


Finally,  the  cell  voltage  can  be  determined  from  the  expression 
as  follows: 


Vcell^O— ^,a  (30) 

where  Vo  is  the  thermodynamic  equilibrium  voltage  of  a  DMFC;  (p®  a 
represents  the  virtual  local  potential  on  the  ribs  of  the  anode  flow- 
field  plate.  Note  that  the  electric  resistance  in  the  electrolyte  phase 
is  implicitly  accounted  by  the  potential  drop  in  the  electrolyte- 
phase-region  when  solving  the  potential  governing  equations. 

Up  to  this  point,  all  the  governing  equations  related  to  the  two- 
phase  mass  transport  and  electrochemical  reactions  have  been 
presented.  In  order  to  exhibit  the  detailed  information  of  the 
governing  equations,  the  computational  sub-domains  for  solving 
variables  are  shown  in  Fig.  3.  The  expressions  of  source  terms  of 
the  specific  equations  and  correlations  of  several  coefficients  are 
listed  in  Table  2. 


(I)  Boundaries  1  and  3:  These  two  boundaries  are  the  interfaces 
between  the  ADL  and  the  ribs  of  the  anode  flow-field  plate. 
These  interfaces  are  impermeable  for  reactants  but  perme¬ 
able  for  electrons: 


dps 

#  =0, 

dx 


dCwv  _  n 
dx 


<Ps  =  <Ps,  a 


(31) 


(II)  Boundary  2:  This  boundary  is  the  interface  between  the  ADL 
and  the  AFC,  which  is  the  inlet  of  methanol  solution  and 
outlet  of  carbon  dioxide  and  water  vapor.  This  boundary  is 
impermeable  for  the  electrons: 


„  „in  r,  „in  ,  ^channel 

Pi  =  Pi, a’  Ps=Pl,a+Pc,a 


s  =  S' 


.channel 


cM  =  cs 


d(ps 


Cwv  =  ^wv’  ^=0 


(32) 


(III)  Boundaries  4  and  8:  The  symmetrical  conditions  for  all 
variables  are  specified  at  these  two  boundaries  as  the  com¬ 
putational  domain  is  a  periodic  unit  of  the  entire  cell: 


dCwv  0  dQ)2  0  9gs  0 
9y  dy  ’  9y 


-f—  =  0  (33) 

ay 


(IV)  Boundaries  5  and  7:  These  two  boundaries  are  the  interfaces 
between  the  DL  and  ribs  at  the  cathode  side.  Similar  to  bound¬ 
aries  1  and  3,  the  conditions  at  these  two  boundaries  can  be 
given  as: 


ac02 

dx 


=  0 , 


dCyw  _ 
dx  ~  ' 


(ps=0 


(34) 


(V)  Boundary  6:  Similar  to  boundary  2,  this  boundary  represents 
the  inlet  of  reactants  on  the  cathode  side.  The  following 
boundary  conditions  at  this  interface  are  specified: 


„  .Jn  ,  ^.channel 
Pl=Pg,c+Pc 


Pg=Pg, 


g,C’ 


^  _  ^channel 


C02  =C“,Cwv  =  0, 


(35) 


(VI)  Boundaries  9  and  12:  As  the  left  and  right  boundaries  of  the 
electrolyte-phase-region,  which  are  impermeable  walls  for 
protons.  Accordingly,  the  flux  of  protons  is  zero  at  these  two 
boundaries: 


Nh+  =  0  (36) 

where  NH+  represents  the  molar  flux  of  protons. 

(VII)  Boundary  10:  This  boundary  is  the  interface  between  the  ACL 
and  MEM.  As  the  membrane  is  treated  as  an  impermeable 
wall  for  the  electrons  and  gas  phase,  the  conditions  at  this 
interface  can  be  given  as: 

Ne-  =  0,  Nwv  =  0  (37) 

(VIII)  Boundary  11:  Similar  to  boundary  10,  conditions  for  this  inter¬ 
face  can  be  specified  as: 


2.7.  Boundary  conditions 


Ne-  =  0,  Nwv  —  0?  Nq2  =  0 


(38) 


Boundary  conditions  are  necessary  for  each  governing  equa¬ 
tion  restricted  to  the  corresponding  computational  sub-domains. 
Conditions  at  different  boundaries  marked  with  Arabic  number 
symbols  are  given  as  follows: 


The  liquid  saturation  at  boundary  2,  s£hannel,  is  one  of  the  key 
parameters  in  the  present  model.  s^hannel  depends  strongly  on  both 
mass  transport  through  the  AFC/ADL  interface  and  the  compli¬ 
cated  two-phase  flow  in  the  AFC.  It  is  a  rather  challenging  work  to 
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Table  2 

Expressions  of  source  terms  and  coefficients  of  the  governing  equations 


Parameters,  symbols 

Expressions  in  sub-domains  of  computational 

ADL 

ACL 

CCL 

CDL 

Generation  rate  of  protons,  im 

0 

ia 

ic-ip 

0 

Generation  rate  of  electrons,  is 

0 

—la 

— (*c — tp ) 

0 

Generation  rate  of  mass  in  gas  phase,  mg 

MWRW 

MwRw  +MCo2Rco2 

Mw^w  +MCo2FCo2  +M02Ro2 

MwRw 

Generation  rate  of  mass  in  liquid  phase,  rhj 

— 

-MWFW  -  MCo2Rqo2 

-MwRw  ~  MCo2Rco2  ~  Mo2Rq2 

—MwRw 

Mole  generation  rate  of  species 

0 

-ia/(6F) 

-ip/(6F) 

0 

Ro2 

0 

0 

— ic/(4F) 

0 

Rwv 

Rw 

Rw 

Rw 

Rw 

Rco2 

0 

ia/(6F) 

iP/(6F) 

0 

Interfacial  molar  transfer  rate  of  water 

kemt(Psw  -yvwPg)^ 

between  liquid  and  vapor  [38],  Fw 

i  b  £(l-s)ywv  f^sat 
+ Kc  2RT  ^W 

-ywvPgXi  -<j) 

Switch  factor 

q=(  1+  |P™~]'wvPs|) 

\  |Pw-ywvps|/ 

Effective  diffusion  coefficients  of  methanol 

i4f  =  DM,|(ss)’'5 

Effective  diffusion  coefficients  of  species  in  gas 

neff_n  T  ■'( 1  -s)-so  1  2 

1  -Ul’Sl  l-e0  J 

phase 

model  the  mass  transport  in  liquid-gas  two-phase  flow  in  the  AFC. 
Additionally,  the  AFC  is  perpendicular  to  and  not  involved  in  the 
computational  domain,  so  some  simplifications  are  necessary.  In 
this  situation,  many  researchers  prefer  to  specify  the  liquid  satura¬ 
tion  in  the  AFC  as  a  constant  [  17,18  ].  In  order  to  reflect  the  saturation 
level  in  the  AFC  for  various  cell  current  density,  the  average  of  liquid 
saturation  at  the  inlet  and  outlet  of  the  AFC  is  applied  in  the  present 
model  instead  of  a  specified  value.  The  mass  conservation  equation 
in  the  AFC  can  be  given  as: 


Din 

lfoutletcoutlet 

T  a 


Vi,afc 


(39) 


uoutlet(1  _s°utlet 


=  V, 


g.afc 


(40) 


where  u°^tlet,  u°uatlet  and  s°utlet  represent  the  liquid  velocity,  gas 
velocity  and  liquid  saturation  at  the  outlet  of  the  AFC.  Qff  is  the 
volumetric  flow  rate  of  methanol  solution  at  the  anode  inlet.  A  is 
the  area  of  the  cross-section  of  the  AFC.  afc  and  Vg  afc  stand  for 
volumetric  flow  rates  of  liquid  and  gas  phases  at  the  AFC/ADL  inter¬ 
face.  According  to  the  present  visualization  research  [21  ],  the  slug 
flow  rather  than  the  bubble  flow  is  more  frequently  encountered  in 
the  AFC.  In  the  slug  flow,  the  large  difference  between  gas  velocity 
and  liquid  velocity  has  significant  influences  on  the  mass  trans¬ 
port  in  the  AFC.  For  more  realistic  modeling  of  the  AFC,  instead 
of  the  homogeneous  model,  a  simplified  expression  derived  from 
the  drift-flux  models  [35,36]  is  applied  to  describe  the  relationship 
between  the  gas  velocity  and  liquid  velocity: 

u°uatlet  =  1 .35  [(1  -  s°utlet)u°uatlet  +  s°utletu°“tlet]  (41 ) 


By  solving  Eqs.  (39)-(41 ),  we  can  obtain  sautlet.  And  then  sahannel 
can  be  calculated  from  the  expression  below: 


cchannel 

^a 


1  +Sgutlet 


2 


(42) 


3.  Numerical  results  and  discussion 

A  self- written  code  based  on  the  finite-volume-method  is  devel¬ 
oped  to  solve  the  governing  equations  iteratively  under  the  baseline 
conditions  and  parameters  given  in  Tables  3  and  4.  Model  validation 
and  the  distributions  of  several  variable  fields  are  present  below. 


3.1.  Model  validation 

Cell  performance  predicted  by  the  present  model  for  the  DMFC 
fed  with  0.25  M  and  0.5  M  methanol  solution  is  compared  with  the 
experimental  data  [37]  in  Fig.  4  under  the  same  operation  condi¬ 
tions,  including  the  anode  flow  rate  of  1.0  ml  min-1  for  15  channels, 
the  cathode  flow  rate  of  1000  ml  min-1  and  the  cell  temperature  at 
75  °C.  In  order  to  validate  the  reliability  of  the  present  model,  the 
same  set  of  parameters  are  used  to  predict  the  polarization  behav¬ 
ior  for  different  methanol  feeding  concentrations.  As  can  be  seen 
from  Fig.  4,  the  predicted  polarization  curves  are  generally  in  good 
agreement  with  the  experimental  data.  The  deviations  of  the  pre¬ 
dicted  data  from  the  experimental  data  in  the  low  current  density 
region  are  probably  caused  by  the  simplifications  used  in  the  mod¬ 
eling  of  the  two-phase  flow  in  the  AFC.  It  can  also  be  found  from  the 
power  density  curves  that  the  highest  power  density  is  achieved  at 
the  front-end  of  the  concentration  polarization  region,  so  it  is  a  wise 
choice  to  let  the  DMFC  work  in  the  Ohm  polarization  region. 

3.2.  Mechanism  analysis  of  species  transport 

The  transport  processes  of  species  in  a  DMFC  are  very 
complicated  as  they  are  intrinsically  coupled  with  both  the  elec- 

Table  3 

Geometric  and  operating  parameters  of  the  DMFC 


Parameters 

Symbols 

Value 

Unit 

Structure  parameters 

Porosity,  thickness 

ADL 

eadl>  ^adl 

0.64,  2.6  x  10-4 

-,  m 

ACL 

eacl>  ^acl 

0.3,  0.2  x  10“4 

-,  m 

MEM 

£mem.  Imem 

0.3, 1.3  x  10“4 

-,  m 

CCL 

£ccl>  ^ccl 

0.3,  0.2  x  10-4 

-,  m 

CDL 

ecdl>  ^cdl 

0.64,  2.6  x  10-4 

-,  m 

Height  of  a  half  rib 

hr 

0.5  x  10-3 

m 

Height  of  channels 

hc 

1  x  10-3 

m 

Length  of  channels 

Lc 

3  x  10-2 

m 

Operation  conditions 

Operation  temperature 

T 

348.15 

I< 

Anode  channel  inlet  pressure 

<a 

101325 

Pa 

Cathode  channel  inlet  pressure 

pin 

101325 

Pa 

Anode  flow  rate 

Qa,in 

1 

ml  min-1 

Inlet  methanol  concentration  at  anode 

250-500 

mol  m3 

Inlet  oxygen  concentration  at  cathode 

cin 

o2 

35 

mol  m3 
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Table  4 

Physicochemical  properties  and  electrochemical  kinetics  parameters  used  in  simulation 
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trochemical  reactions  and  the  mass  transport  between  liquid  and 
gas  phases.  So  better  understanding  of  fundamental  species  trans¬ 
port  mechanisms  in  a  DMFC  is  essential  to  the  optimization  of  the 
cell  design. 


Fig.  4.  Comparison  of  cell  performance  predicted  by  the  present  model  and  exper¬ 
imental  data  [37]. 


The  concentration  distributions  of  methanol  and  oxygen  are 
given  in  Fig.  5(a)  and  (b)  respectively.  A  sharp  decrease  in  methanol 
concentration  can  be  seen  from  the  AFC  to  the  ACL  due  to  the  con¬ 
sumption  in  the  MOR.  The  methanol  concentration  also  decreases 
from  the  under-channel  region  to  the  under-rib  region  because  ribs 
of  the  flow-field  plate  limit  the  access  of  methanol  solution  to  the 
ACL.  As  a  contrast,  the  distribution  of  oxygen  concentration  in  the 
CDL  and  CCL  is  nearly  uniform,  which  can  be  observed  from  Fig.  5(b). 
The  largest  relative  difference  of  oxygen  concentration  is  0.63%.  This 
is  because  the  diffusion  coefficient  of  oxygen  is  several  orders  of 
magnitude  higher  than  that  of  methanol  and  the  density  of  the  gas 
phase  is  very  low.  Fig.  5  indicates  that  the  degradation  of  cell  perfor¬ 
mance  caused  by  the  concentration  polarization  on  the  electrodes 
is  mainly  due  to  the  mass  transport  limitation  of  methanol.  In  sec¬ 
tion  3.4,  the  non-uniform  distribution  of  MOR  rate  mainly  caused 
by  the  non-uniform  distribution  of  methanol  concentration  in  the 
ACL  will  also  be  observed. 

The  total  flux  of  methanol  transport  consists  of  two  contribu¬ 
tions  with  different  transport  mechanisms:  diffusion  caused  by  the 
gradient  of  methanol  concentration  as  can  be  seen  from  Fig.  5(a) 
and  convection  caused  by  the  bulk  motion  of  the  liquid  phase.  The 
liquid  motion  in  the  MEA  is  driven  by  the  gradient  of  liquid  pres¬ 
sure  and  the  electro-osmotic  drag  force.  Noteworthy  is  that  most 
of  the  previous  models  for  DMFCs  mainly  focus  on  the  transport 
phenomena  in  the  anodic  or  cathodic  porous  media.  The  present 
model  also  provides  a  modeling  of  the  simultaneous  methanol  and 
water  transport  across  the  MEM.  Fig.  6  shows  the  distributions  of 
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Fig.  5.  Distribution  of  (a)  methanol  concentration  in  four  layers:  ADL,  ACL,  MEM  and  CCL,  (b)  oxygen  concentration  in  the  cathode  (Cell  voltage:  0.11  V). 


velocity  fields  of  liquid  and  gas  phases  in  the  upper  half  of  the  com¬ 
putational  domain.  Apparently,  the  complex  liquid-gas  two-phase 
counter-flow  behaviour  can  be  seen  in  the  porous  MEA.  As  shown 
in  Fig.  6(a),  the  methanol  solution  moves  from  the  AFC  to  the  ACL  in 
two  directions:  directly  to  the  under-channel  region  of  the  ACL  and 
snaking  through  the  ADL  to  the  under-rib  region  of  the  ACL.  And 
then,  the  excess  water  and  methanol  in  the  ACL  permeate  through 
the  MEM  to  the  CCL,  where  methanol  is  totally  oxidized  and  the 
ORR  takes  place  to  form  water.  The  produced  water,  along  with 
the  water  crossing  over  from  the  anode,  transfers  through  the  CDL 


0  0.00069 


X(m) 


to  the  CFC.  By  contrast,  the  gas  phase  movement  in  the  opposite 
direction  can  be  seen  from  Fig.  6(b),  as  oxygen  is  consumed  in  the 
CCL  while  carbon  dioxide  is  simultaneously  generated  in  the  ACL. 
Because  the  MEM  is  regarded  as  a  gas  insulator,  there  is  a  blank  of 
the  gas  velocity  field  in  this  region. 

In  order  to  quantify  the  effects  of  diffusion  and  convection  on  the 
methanol  transport  in  the  ACL,  their  contributions  to  the  mean  cur¬ 
rent  density  for  0.25  M  and  0.5  M  methanol  feeding  concentrations 
are  compared  in  Fig.  7.  It  can  be  clearly  seen  that  convection  only 
accounts  for  a  small  fraction  of  the  total  flux  of  methanol.  More  than 


(b)  5E-3  m  s'1 


Fig.  6.  Distribution  of  velocity  vector  of  (a)  liquid  phase  and  (b)  gas  phase  in  the  upper  half  of  the  computational  domain  (Cell  voltage:  0.11  V). 
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Cell  current  density  (A  rrf2) 


Fig.  7.  Share  of  contributions  of  diffusive  and  convective  transport  to  the  cell  current 
density  for  various  cell  voltage. 


Fig.  8.  Liquid  pressure  along  x  direction  for  different  locations  on  y  direction  (Cell 
voltage:  0.11  V). 


90%  of  the  mean  current  density  results  from  methanol  diffusive 
transport,  which  indicates  that  the  mass  transport  of  methanol  in 
the  ACL  is  dominated  by  diffusion  mechanism.  Besides  that,  another 
trend  can  be  observed  that  the  percentage  of  convection  for  a  given 
methanol  feeding  concentration,  0.25  M  or  0.5  M,  decreases  with 
an  increase  in  cell  current  density.  This  phenomena  result  from 
the  combined  effects  of  liquid  velocity  and  methanol  concentration 
on  methanol  transport.  When  the  current  density  becomes  larger, 
more  methanol  and  water  will  be  consumed  in  the  ACL,  and  this 
will  lead  to  an  increase  in  the  liquid  velocity  and  a  decrease  in  the 
methanol  concentration.  Accordingly,  it  is  easy  to  conceive  that  the 
contribution  of  the  increased  liquid  velocity  to  the  methanol  con¬ 
vection  is  weakened  or  even  submerged  by  the  decreased  methanol 
concentration  while  the  diffusive  transport  contributes  more  sig¬ 
nificantly  to  the  total  methanol  flux  due  to  the  increase  in  the 
concentration  gradient  of  methanol.  According  to  the  above  anal¬ 
ysis,  it  can  be  seen  that  convection  is  a  weak  point  in  methanol 
transport  processes.  Efforts  on  enhancing  the  convective  methanol 
transport  can  significantly  improve  the  cell  performance. 

3.3.  Analysis  of  water  and  methanol  crossovers  through  the  MEM 

Water  and  methanol  crossovers  through  the  MEM  from  the 
anode  to  the  cathode  are  two  of  the  key  technological  challenges 
in  the  research  of  DMFCs.  The  total  flux  of  methanol  crossover  is 
comprised  of  three  contributions,  namely:  diffusion  due  to  the  con¬ 
centration  gradient  of  methanol,  back  convection  due  to  the  bulk 
motion  under  adverse  pressure  gradient  and  electro-osmotic  drag 
flow  due  to  the  transport  of  protons  in  the  electrolyte  phase.  As  it 
is  assumed  that  the  MEM  is  fully  hydrated  and  no  water  diffusion 
occurs,  water  crossover  mainly  rests  with  the  back  convection  and 
electro-osmotic  drag  flow. 

Fig.  8  shows  the  liquid  pressure  profiles  along  x  direction  of  the 
computational  domain  at  three  different  locations  on  y  direction.  It 
can  be  seen  that  the  liquid  pressure  increases  rapidly  in  the  MEM. 
This  huge  adverse  pressure  gradient  causes  the  back  convection  of 
methanol  solution  from  the  cathode  to  the  anode,  which  can  reduce 
methanol  crossover,  and  is  beneficial  to  the  cell  performance.  In 
this  model,  we  considered  the  contribution  of  electro-osmotic  drag 
to  liquid  velocity  not  only  in  the  MEM  but  also  in  CLs.  It  can  be 
observed  from  the  detailed  view  of  the  liquid  pressure  profiles  that 
the  driving  force  of  liquid  velocity  is  transferred  gradually  from 
pressure  gradient  to  electro-osmotic  drag  force,  exhibited  as  the 


minimum  and  maximum  points  of  the  liquid  pressure  profiles.  This 
proves  that  it  is  necessary  to  consider  the  electro-osmotic  drag  flow 
in  CLs  in  order  to  avoid  the  overestimate  of  the  contribution  of 
pressure  gradient  to  the  liquid  velocity,  especially  on  the  cathode 
side.  Also  in  Fig.  8,  the  liquid  pressure  in  the  under-rib  region  of 
the  CDL  is  higher  than  that  in  the  under-channel  region.  This  trend 
agrees  well  with  the  results  reported  in  [15,17]  as  the  liquid  water 
is  prone  to  accumulate  in  the  corner  under  ribs.  However,  a  higher 
and  sharply  decreasing  liquid  pressure  in  the  under-channel  region 
in  the  CCL  can  also  be  seen  from  the  detailed  view  of  liquid  pressure 
in  the  CCL.  This  pressure  profile  indicates  that  the  water  flux  in  the 
under-channel  region  of  the  CCL  is  heavier  than  that  in  the  under¬ 
rib  region.  The  water  flux  in  the  CCL  includes  water  generated  by 


Fig.  9.  Distribution  of  potential  of  electrolyte  phase  and  the  local  current  of  protons 
in  the  electrolyte-phase-region  (Cell  voltage:  0.11  V). 
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Cell  voltage  (V) 

Fig.  10.  Percentage  of  contributions  resulting  from  different  methanol  transport 
mechanisms  to  the  total  methanol  crossover  flux  for  various  cell  voltage. 

the  ORR  and  more  importantly,  flux  of  water  crossover  through  the 
MEM,  which  is  the  larger  part  of  the  total  water  amount. 

The  distributions  of  potential  and  the  local  current  of  protons  in 
electrolyte-phase-region  are  shown  in  Fig.  9.  Similar  to  the  distri¬ 
bution  of  methanol  concentration,  the  electrolyte  phase  potential 
decreases  from  the  ACL  to  the  CCL  and  from  the  under-channel 
region  to  the  under-rib  region.  This  is  a  result  of  the  higher  proton 
generation  rate  in  the  under-channel  region  of  the  ACL  as  protons 
move  from  the  ACL  to  the  CCL  driven  by  the  electrical  field  force.  The 
generation  and  consumption  of  protons  in  terms  of  the  increase  and 
decrease  in  the  local  current  vector  can  also  be  observed  from  the 
sub-regions:  x  =  0.00026-0.00028  m  and  x  =  0.00041  -0.00043  m. 

Figs.  8  and  9  provide  a  detailed  description  of  water  crossover 
through  the  MEM.  A  common  feature  of  these  two  figures  is  that 
both  the  liquid  pressure  gradient  and  electro-osmotic  drag  force 
promote  the  migration  of  protons,  methanol  and  water  from  the 
under-channel  region  to  the  under-rib  region  in  the  electrolyte- 
phase-region,  which  can  also  be  clearly  seen  in  Fig.  6(a).  This 
migration  process  can  reduce  the  non-uniformity  of  the  DMFC, 
make  better  usage  of  the  catalyst  in  the  CCL  and  thus  is  beneficial 
to  the  cell  performance. 

In  order  to  graphically  show  the  effects  of  each  methanol  trans¬ 
port  mechanism  on  the  methanol  crossover,  the  percentage  of  their 
contributions  to  the  total  flux  of  methanol  crossover  are  calcu¬ 
lated.  As  can  be  seen  in  Fig.  10,  the  scale  of  the  percentage  is  set 
to  begin  with  60%  to  give  prominence  to  the  difference  between 
cases  for  various  cell  voltages.  Clearly,  the  diffusive  transport  dom¬ 
inates  methanol  crossover  with  a  share  of  more  than  85%.  The  back 
convection  of  methanol  solution  can  reduce  methanol  crossover 
flux  by  approximately  2%.  So  carefully  adjusting  the  inlet  pressure 
of  the  AFC  and  CFC  or  making  some  structure  improvements,  such 
as  adding  a  micro-pore  layer  between  the  CCL  and  CDL  will  sup¬ 
press  methanol  crossover  by  increasing  the  back  convection  flux  of 
methanol  solution. 

3.4.  Analysis  of  the  distributions  of  electrochemical  reaction  rates 
in  CLs 

Electrochemical  reactions  occurring  in  CLs  of  the  DMFC  are 
very  complicated  multi-step  processes  which  are  always  described 
by  empirical  expressions  like  Eqs.  (19)  and  (20).  Because  the 
overpotential  is  at  the  power  position  of  the  Euler’s  constant 
(approximately  2.7183),  it  has  significant  impact  on  the  electro¬ 
chemical  reaction  rates.  In  the  present  model,  we  calculated  the 


Fig.  11.  Distribution  of  overpotential,  MOR  rate  and  methanol  concentration  in  the 
ACL  (Cell  voltage:  0.31  V). 

distribution  of  overpotential  and  the  effect  of  overpotential  on  elec¬ 
trochemical  reaction  rates,  as  shown  in  Figs.  11-13.  It  should  be 
noted  that  the  non-dimensional  relative  scale  is  adopted  for  the 
investigated  variables  to  clearly  show  their  distributions  because 
the  methanol  concentration  and  electrochemical  reaction  rates 
vary  within  several  orders  of  magnitude.  Values  of  variables  at  point 
(x  =  0.00026  m,  y  =  0.001  m)  are  set  as  the  reference  values  in  the 
ACL,  and  values  of  variables  at  point  (x  =  0.00041  m,  y  =  0.001  m) 
are  set  as  the  reference  values  in  the  CCL. 


Fig.  12.  Distribution  of  oxygen  concentration,  overpotential,  ORR  rate  and  PCR  rate 
in  the  CCL  (Cell  voltage:  0.31  V). 


Z.  Miao  et  al.  /  Journal  of  Power  Sources  185  (2008)  1233-1246 


1243 


The  distributions  of  anode  overpotential,  MOR  rate  and 
methanol  concentration  in  the  ACL  are  presented  in  Fig.  11. 
Apparently,  methanol  concentration  decreases  sharply  from  the 
under-channel  region  to  the  under-rib  region  and  gently  along  x 
direction,  while  the  distribution  of  anode  overpotential  displays  a 
reverse  but  more  gradual  trend.  As  a  result  of  the  combined  effects 
of  methanol  concentration  and  anode  overpotential,  the  distribu¬ 
tion  of  MOR  rate  appears  in  the  shape  of  a  saddle,  similar  to  the 
distribution  of  methanol  concentration  but  having  an  uprising  in 
the  region  close  to  the  MEM.  So  it  can  be  concluded  that  the  MOR 
rate  in  the  ACL  is  mainly  dominated  by  the  methanol  concentra¬ 
tion  while  the  work  point  of  the  cell  is  within  the  mass  transport 
controlled  region.  Also  in  Fig.  11,  the  anode  overpotential  exhibits 
a  self-regulation  characterize  to  the  effect  of  non-uniform  distri¬ 
bution  of  methanol  concentration  on  the  MOR  rate  by  making  the 
distribution  of  MOR  rate  more  uniform  throughout  the  ACL. 

It  has  been  mentioned  about  the  results  in  Fig.  5  that  the  dis¬ 
tribution  of  oxygen  concentration  in  the  CCL  is  almost  uniform, 
which  can  also  be  seen  in  Fig.  12.  If  the  overpotential  in  the  CCL 
is  set  as  a  constant,  the  uniform  distribution  of  ORR  rate  will  be 
expected  logically.  Flowever,  an  extremely  non-uniform  distribu¬ 
tion  of  ORR  rate  can  be  seen  in  Fig.  12  due  to  the  slightly  decrease 
in  the  cathode  overpotential  in  the  under-rib  region.  These  distri¬ 
butions  indicate  that  the  non-uniformity  of  the  ORR  rate  in  the  CCL 
depends  highly  on  the  cathode  overpotential.  As  mentioned  above, 
the  total  amount  of  oxygen  consumed  in  the  ORR  includes  two 
parts:  oxygen  reacted  with  protons  transported  from  the  anode, 
named  proton  consumption  reaction  (PCR)  in  this  paper,  and  oxy¬ 
gen  consumed  in  the  MDOR.  In  Fig.  12,  a  large  gap  between  the  PCR 
rate  and  the  ORR  rate  in  the  region  near  the  MEM  can  be  observed, 
and  this  indicates  that  although  most  of  the  oxygen  in  the  CCL  is 
consumed  in  the  region  near  the  MEM,  only  part  of  the  oxygen  par¬ 
ticipates  in  the  PCR.  And  this  situation  will  get  worse  when  the 
DMFC  works  at  a  lower  current  density  as  more  methanol  will  per¬ 
meate  through  the  MEM.  On  the  other  hand,  Fig.  12  also  predicts 
a  non-uniform  distribution  of  the  MDOR  rate  which  can  be  clearly 


Fig.  13.  Distribution  of  overpotential,  methanol  concentration  and  MDOR  rate  in 
the  CCL  (Cell  voltage:  0.31  V). 


seen  in  Fig.  13.  The  distribution  of  methanol  concentration  in  the 
CCL  is  extremely  non-uniform  and  the  MDOR  rate  decreases  more 
sharply  than  methanol  concentration  along  both  x  and  y  directions 
because  of  the  effect  of  the  slight  non-uniform  distribution  of  cath¬ 
ode  overpotential.  The  results  showed  in  Figs.  11-13  indicate  that 
accounting  for  the  simultaneous  transport  processes  of  protons, 
methanol  and  water  in  the  electrolyte-phase-region  is  helpful  in 
achieving  a  more  realistic  prediction  of  the  distribution  of  electro¬ 
chemical  reaction  rates  in  CLs. 

3.5.  Effect  of  the  agglomerate  radius  and  overlapping  angle 

A  modified  agglomerate  model  is  introduced  in  this  work  to 
consider  the  effect  of  microstructure  of  the  CLs  on  cell  perfor¬ 
mance.  Fig.  14  presents  the  effect  of  the  agglomerate  radius  on 
polarization  curves.  The  Ragg  =  0  case  corresponds  to  the  homoge¬ 
neous  model  in  which  the  transport  resistance  in  the  agglomerates 
is  not  taken  into  account.  It  can  be  seen  that  the  homogeneous 
model  overestimates  the  cell  performance  due  to  neglect  of  the 
mass  transfer  resistance  induced  by  agglomerates.  Fig.  14  also 
shows  that  the  cell  voltage  decreases  with  an  increase  in  the 
agglomerate  radius  because  of  the  higher  mass  transfer  resistance 


Fig.  15.  Effect  of  the  overlapping  angle  of  agglomerates  on  the  cell  performance. 
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in  larger  agglomerate.  Results  in  Fig.  14  indicate  that  fully  dis¬ 
persion  of  catalyst  particles  in  the  fabrication  processes  of  the 
electrodes  to  achieve  smaller  size  can  lead  to  a  better  cell  perfor¬ 
mance. 

Besides  that,  the  overlapping  angle,  another  key  parameter  of 
the  agglomerate,  can  reflect  the  topological  of  the  agglomerates 
to  a  certain  extent.  The  effect  of  the  overlapping  angle  on  cell 
performance  is  shown  in  Fig.  15.  Note  that  agglomerates  with 
overlapping  angles  equal  to  0  and  90  degrees  correspond  to  the 
spherical  and  cylindrical  agglomerate  geometry  respectively. 
And  our  modified  model  can  predict  the  general  microstructure 
of  agglomerate  with  an  overlapping  angle  between  0  and  60 
degrees.  An  increase  in  the  overlapping  angle  only  leads  to  a  slight 
decrease  in  the  cell  voltage,  implying  that  the  overlapping  angle 
in  the  present  agglomerate  model  has  a  rather  little  effect  on  the 
polarization  curve.  This  little  effect  is  probably  attributed  to  the 
decrease  in  the  specific  surface  area  of  the  agglomerates  when  the 
overlapping  angle  becomes  larger. 


4.  Conclusions 
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Appendix  A.  Derivation  of  the  correction  factor  §o2  based 
on  the  modified  agglomerate  model 


It  can  be  seen  from  Fig.  2  that  in  order  to  supply  access 
for  the  electrons,  the  neighboring  agglomerates  should  share  an 
overlapped  region  with  each  other.  So  a  segmental  spherical 
agglomerate  geometry  is  adopted  in  this  model  with  a  simplifi¬ 
cation  that  only  the  transfer  resistance  on  the  radial  direction  is 
considered.  The  conservation  equation  of  oxygen  concentration  in 
sub-domain  I:  the  Nation  coating  region,  can  be  given  as 


Ragg  COS  0 
r  +  dr 


4jt{r  +  dr)2DN 


d_ 

dr 


dC o2,n 
dr 


A  two-dimensional,  two-phase  mass  transport  model  adopting 
a  modified  agglomerate  model  has  been  developed  to  investigate 
the  mechanisms  of  species  transport  processes  and  the  distribu¬ 
tions  of  several  key  variables  in  the  DMFC.  The  microstructure  of 
the  CLs  and  the  effect  of  electro-osmotic  drag  on  the  liquid  motion 
in  the  electrolyte-phase-region  are  considered.  The  electrons  and 
protons  transport  processes  are  also  taken  into  account  to  numer¬ 
ically  calculate  the  distribution  of  overpotential  in  both  ACL  and 
CCL.  Based  on  the  analysis  of  the  numerical  results,  conclusions  are 
summarized  as  follows: 

(1 )  The  distribution  of  methanol  concentration  is  significantly  non- 
uniform  along  both  x  and  y  directions  in  the  computational 
domain.  As  a  contrast,  almost  uniform  distribution  of  oxy¬ 
gen  concentration  appears  in  the  cathodic  porous  region  due 
to  the  lower  transport  resistance  for  gas  phase  in  the  CDL. 
Compared  with  convection,  diffusive  transport  dominates  the 
methanol  transport  and  has  major  effect  on  the  working  current 
density. 

(2)  Methanol  crossover  mainly  results  from  the  diffusive  transport 
while  water  crossover  mainly  rests  with  the  electro-osmotic 
drag  flow.  In  the  electrolyte-phase-region,  methanol  and  water 
is  prone  to  transfer  from  the  under-channel  region  to  the 
under-rib  region  in  addition  to  the  main  transport  processes 
along  x  direction.  Numerical  results  also  show  that  it  is  nec¬ 
essary  to  consider  the  migration  process  of  protons  in  the 
electrolyte  phase  to  achieve  a  detailed  description  of  the  trans¬ 
port  processes  of  protons,  methanol  and  water  throughout  the 
electrolyte-phase-region. 

(3)  The  distributions  of  different  variables  interact  with  each  other 
and  this  convoluted  relationship  can  be  mainly  understood 
as  that  non-uniform  distribution  of  methanol  concentration 
results  in  non-uniform  distribution  of  MOR  rate,  and  the  lat¬ 
ter  will  lead  to  a  non-uniform  distribution  of  overpotential  in 
the  CCL.  An  extremely  non-uniform  distribution  of  ORR  rate 
appears  in  the  CCL  because  of  the  significant  effect  of  cathode 
overpotential  on  the  ORR  rate  even  though  the  distribution  of 
oxygen  is  nearly  uniform. 

(4)  The  numerical  results  indicate  that  the  modified  agglomerate 
model  in  this  paper  can  reflect  a  more  reasonable  and  realistic 
microstructure  of  the  CCL  as  well  as  the  effect  of  this  microstruc¬ 
ture  on  the  cell  performance.  The  results  also  suggest  that  a 
smaller  agglomerate  size  is  beneficial  to  the  cell  performance. 
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=  0 
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By  simplifying  Eq.  (A.l ),  Eq.  (23)  can  be  obtained 

d2C0  N+1d^N=(] 
dr2  r  dr 

In  sub-domain  II,  the  conservation  equation  can  be  written  in 
the  form  similar  to  Eq.  (23)  with  a  non-zero  source  term: 


d2Co2,agg  ,  1  dQ)2,agg  _  kv  r 

dr2  r  dr  ~  Deff  °2,agg 
02,agg 


(24) 


The  geometry  of  sub-domain  III  is  the  same  as  the  spheri¬ 
cal  agglomerate,  so  Eq.  (25)  with  the  same  form  as  the  standard 
concentration  conservation  equation  in  a  spherical  coordinate  is 
suitable  to  formulate  the  mass  transfer  process  in  this  region: 


d2Co2,agg  (  2  dCo2,agg  _  kv  r 

dr2  r  dr  ~  D^ff  °2,agg 
02,agg 


(25) 


Three  variables:  Co2,i,  Co2,n  and  Co2,m  are  defined  in  this  model 
to  represent  the  oxygen  concentration  at  the  outer  surfaces  of  sub- 
domains  I— III.  According  to  the  spherical  agglomerate  model,  we 
can  obtain: 
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where  the  Thiele  modulus  is  defined  as 
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Eqs.  (A.4)  and  (A.5)  is  used  as  the  boundary  condition  of  Eq. 
(24),  which  is  a  Bessel  function  and  has  a  power  series  solution  as 
follows: 


Co2,agg(r)  =  Co2,i^Q(r  -  Rin)n  Rin  <  r  <  Ragg  (A.6) 

n=0 
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where  the  coefficients  can  be  given  as 

Co  =  i,  c;  =  2_(0S  coth  4>s  - 1) 

*Mn 

(A.7) 

r  Yfk=0^k  +  Pan-kCk+3  +  b0Ck 

(n  +  l)(n  +  2) 

(A.8) 

ak=  \  '  k  =  0,1,2,... 

C 

(A.9) 

(A.10) 

The  information  at  the  outer  interface  of  sub-domain  II  which 
can  be  obtained  from  Eq.  (A.6)  is  as  follows: 
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n= 0 
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n= 1 


The  correction  factor  §agg  in  view  of  the  transfer  resistance  in 
the  agglomerate  is  defined  as 


5aggE)o2 ,  agg  ( ^°2 ,  agg  /^r  )  |  r= 

VaggkvCo2,u 


(A.13) 


where  Sagg  and  Vagg  represent  the  surface  area  and  the  volume  of  the 
agglomerate.  By  substituting  Eqs.  (A.ll)  and  (A.12)  into  Eq.  (A.13), 
the  expression  of  §agg  is  rewritten  as 
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W  =  ^[C'(Ragg-Rin)n]  (A.15) 

n= 0 

/((9)  = - — 4 -  0°  <  e  <  60°  (A.16) 

[2  -  (1  -  cos  0)2{ 2  +  cos  6)] 

For  sub-domain  I,  the  solution  of  Eq.  (23)  can  be  expressed  as: 


h  height  of  channel  or  rib  (m) 

i  electrochemical  reaction  rate  (A  m-3 ) 

I  current  density  (A  m-2 ) 

I  current  vector  (A  m-2 ) 

kc  condensation  rate  of  water  (s-1 ) 

ke  evaporation  rate  of  water  (atm-1  s-1 ) 

/<h  Henry’s  law  constant 

krg  relative  permeability  of  gas  phase 

kr i  relative  permeability  of  liquid  phase 

K  absolute  permeability  of  porous  media  (m2) 

I  length  of  the  channel  (m) 

m  source  terms  in  mass  conservation  equations  (kg  m-3  s-1 ) 

M  molecular  weight  (kg  mol-1 ) 

nd  electro-osmotic  drag  coefficient 

N  molar  flux  (mol  m2  s-1) 

p  pressure  (Pa) 

pc  capillary  pressure  (Pa) 

q  switch  factor 

Q  volume  flow  rate  (ml  min-1 ) 

R  gas  constant  (J  mol-1 1<-1 ) 

R  source  term  in  species  conservation  equations 

(molm-3  s-1) 

interfacial  transfer  rate  of  water  (mol  m-3  s-1 ) 
s  liquid  saturation 

T  temperature  (K) 

u  superficial  velocity  vector  (m  s-1 ) 

Vo  thermodynamic  equilibrium  voltage  (V) 

Vcell  cell  voltage  (V) 

Greek  symbols 

a  transfer  coefficient 

y  reaction  order 

5n  thickness  of  Nafion  coating  (m) 

e  porosity  of  the  porous  media 

k  conductivity  of  membrane  phase  (£2-1  m-1 ) 

fi  viscosity  (kg  m-1  s-1 ) 

0  overlapping  angle  (°) 

6C  contact  angle  (°) 

p  density  (kg  m-3) 

a  interfacial  tension  (N  m-1 ) 


Co2,N(r)  =  a  In  r  +  b  (A.17) 

where  a  and  b  are  the  coefficients  related  to  the  boundary  con¬ 
ditions  at  the  interface  between  sub-domains  I  and  II  and  can  be 
derived  by  considering  Eqs.  (A.ll )  and  (A.12).  So  the  correction  fac¬ 
tor  £n  in  view  of  the  transfer  resistance  in  the  Nafion  coating  can 
be  given  as 


1  +  (2/3)(D|gg/D]\[)(0s  coth  (ps  —  1)  In  (flagg  +  <$N)/l^agg 

Finally,  the  general  correction  factor  £o2  concerning  the  overall 
mass  transport  processes  of  oxygen  from  the  gas  pore  to  the  active 
sites  of  catalyst  particles  can  be  calculated  by  the  expression  below: 

£o2  =  £^-5  N^agg  (A.19) 

Appendix  B.  Nomenclature 


List  of  symbols 

A  specific  area  in  the  catalyst  layer  (m2  m-3 ) 

C  concentration  ( mol  m-3 ) 

D  diffusivity  (m2  s-1 ) 

F  Faraday  constant  (96485  C  mol-1 ) 


Superscripts 

channel 

the  flow  channel 

eff 

effective  value 

in 

inlet  of  the  flow  channel 

ref 

reference  value 

sat 

saturated 

Subscripts 

a 

anode 

acl 

anode  catalyst  layer 

adl 

anode  diffusion  layer 

agg 

the  agglomerate 

c 

cathode 

ccl 

cathode  catalyst  layer 

cdl 

cathode  diffusion  layer 

e- 

electrons 

g 

gas  phase 

H+ 

protons 

in 

the  inner  sphere 

1 

liquid  phase 

m 

the  membrane  phase 

mem 

membrane 

N 

Nafion  phase 

02 

oxygen 
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s  the  solid  phase 

WV  water  vapor 
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